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ABSTRACT 

We use numerical simulations to study the evolution of triaxial elliptical galaxies with central black 
holes. In contrast to earlier studies which used galaxy models with central density "cores," our galaxies 
have steep central cusps, as observed in real ellipticals. As a black hole grows in these cuspy triaxial 
galaxies, the inner regions become rounder owing to chaos induced in the orbital families which populate 
the model. At larger radii, however, the models maintain their triaxiality, and orbital analyses show 
that centrophilic orbits there resist stochasticity over many dynamical times. While black hole induced 
evolution is strong in the inner regions of these galaxies, and reaches out beyond the nominal "sphere 
of influence" of a black hole, our simulations do not show evidence for a rapid global transformation of 
the host. The triaxiality of observed elliptical galaxies is therefore not inconsistent with the presence of 
supermassive black holes at their centers. 

Subject headings: galaxies: elliptical, galaxies: kinematics and dynamics, galaxies; structure, methods: 
N-body simulations 



L INTRODUCTION 

Observations indicate that most, and perhaps all ellipti- 
cal galaxies harbor supermassive black holes at their cen- 
ters (e.g. Gebhardt et al. 2000; Richstone et al. 1998; 
but see Gebhardt et al. 2GGlb). In fact, best-fit mod- 
els of black hole demography indicate that roughly 97% 
of ellipticals contain such black holes (Magorrian et al. 
1998). The masses of these black holes seem to be cor- 
related with properties of the host bulge; current dy- 
namical estimates have yielded black hole masses of or- 
der O.OOSMbuigc (Kormendy & Richstone 1995; Magor- 
rian et al. 1998; van der Marel 1999). There is also 
apparently a trend between black hole mass and galaxy 
velocity dispersion, implying that there is a Fundamen- 
tal Plane even in the four-dimensional space defined by 
[logMBH, logi, log cTe, logi?e] (Gebhardt ct al. 2001a; Fer- 
rarese & Merritt 2001). These correlations suggest that 
galaxy formation and the formation of central black holes 
are deeply connected. 

It is believed that a central supermassive black hole can 
profoundly influence the evolution of its host galaxy. Mas- 
sive black holes dominate the galactic potential inside a ra- 
dius of influence tbh ~ 100(AfBH,9/o'o 200) P*^ '^li^re Mbh,9 
is the black hole mass in units of lO^M© and cro,20o is the 
central velocity dispersion in units of 200 km/sec. Within 
^ < ''BH, the three-dimensional structure and phase space 
of a galaxy will be determined largely by the black hole. 
Moreover, the effects of a central black hole can reach far 
beyond tbh- Owing to discrete encounters with individ- 
ual stars, central black holes will "wander," in a manner 
akin to molecular Brownian motion (see, e.g. Chatterjee, 
Hernquist & Loeb [2001] for a recent analysis of this phe- 
nomenon), effectively increasing tbh- During galaxy can- 



nibalism, supermassive black holes disrupt even dense low- 
mass companions in radial encounters (HoUey-Bockelmann 
& Richstone 2000; Merritt & Cruz 2001) and may tidally 
torque debris into a nuclear disk following mergers from 
eccentric orbits (HoUey-Bockelmann & Richstone 2000). 
Black hole binaries created during galaxy mergers may ex- 
plain the fiat density profiles seen in the inner regions of 
the largest ellipticals (Makino & Ebisuzaki 1996; Quinlan 
& Hernquist 1997). 

Central black holes can induce secular evolution as well. 
According to one long-standing suggestion, the growth of 
a massive black hole in a triaxial potential can destabilize 
centrophilic box orbits through stochastic diffusion, driv- 
ing the global shape of a galaxy towards axisymmetry in 
a few crossing times (Gerhard & Binney 1985; Norman, 
May, & van Albada 1985; Merritt & Quinlan 1998; Wach- 
hn & Ferraz-Mello 1998; Valluri & Merritt 1998). If this 
picture of black hole induced evolution away from triaxial- 
ity is correct, it will have serious ramifications for our un- 
derstanding of the fueling of active galactic nuclei (AGN) 
and the observed correlations between the properties of 
black holes and their host galaxies. Indeed, such evolu- 
tion could lead to a self-regulation of AGN activity and 
black hole growth in ellipticals. If triaxiality is linked to 
AGN fueling, as a consequence of gas being unable to set- 
tle into closed orbits in such a potential (Norman & Silk 
1983), then an evolution towards axisymmetry would nat- 
urally lead to a reduction in the feeding of a central black 
hole. Such a scenario might also explain, at least in part, 
the correlation between black hole and bulge mass in early 
type galaxies (Valluri & Merritt 1998). 

While theoretical arguments paint a compelling picture 
for a link between central black holes and the loss of tri- 
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axiality in galaxies, observational support remains some- 
what problematic, owing to the difficulty of inferring the 
three dimensional structure of galaxies from their pro- 
jected properties. However, recent studies indicate that at 
least some luminous ellipticals are triaxial (Binney 1976; 
Franx et al. 1991; Tremblay & Merritt 1995; Bak & Statler 
2000). The triaxiality of these systems indicates either 
that they contain at most low mass black holes (insuffi- 
cient to drive evolution) or that the coupling between black 
holes and galaxy shape is not well understood. The for- 
mer possibility appears to be inconsistent with the notion 
that luminous AGN are preferentially found in luminous 
ellipticals. 

These considerations emphasize the need for realistic 
models of galaxies containing central black holes. Earlier 
theoretical studies have employed either triaxial models 
with constant density cores (Norman et al. 1985; Merritt 

6 Quinlan 1998) or spherical models with central density 
cusps (e.g. Sigurdsson et al. 1995). Neither of these two 
choices is entirely satisfactory. Observations of ellipticals 
made with the Hubble Space Telescope show that these 
galaxies typically have central density cusps oc r~'^ with 

7 = 0.5 — 2 (Lauer et al. 1995). Galaxies with density 
cusps support different stellar orbits than e.g. 7 = core 
models (Gerheard & Binney 1985; Gerhard 1986; Pfen- 
niger & de Zeeuw 1989; Schwarzschild 1993; de Zeeuw 
1996; Merritt 1999; Holley-Bockelmann et al. 2001), alter- 
ing the response of a galaxy to the growth of a black hole. 
Likewise, spherical galaxies will support different orbital 
families than triaxial ones, whether or not they have cen- 
tral density cusps. While simulations employing spherical 
models have demonstrated the growth of black hole in- 
duced cusps and the polarization of the velocity ellipsoid 
in such objects (Sigurdsson et al. 1995), they of course 
have little to say about the evolution of triaxial ellipticals. 

To address these shortcomings, we present a study of 
black hole growth in triaxial elliptical galaxies with central 
density c;iisps, using the "adiabatic squeezing" technique 
described in Holley-Bockelmann et al. (2001; hereafter, 
"Paper 1" ) to generate models with prescribed shapes and 
cusps that are stable for many dynamical times. Using N- 
body simulations, we then adiabatically grow a black hole 
in these galaxies over sevciral crossing times. We charac- 
terize the orbital families which populate the models, and 
find that many of the highly bound boxes, boxlets, and ec- 
centric tubes are transformed into chaotic orbits a clear 
signature of the influence of the spherical central potential 
in a larger, nearly unchanged, triaxial figure. 

The paper is organized as follows. Section 2 summa- 
rizes our technique for generating a stable triaxial model 
of a galaxy, with subsequent growth of a central black 
hole. Section 3 shows the outcome of black hole growth 
in a model with central density cusp 7 = 1.0 and initial 
half mass axis ratios b/a, c/a = 0.85, 0.7, and discusses the 
structural changes occurring in this model as a result of the 
growth of a black hole of mass Mbh = O.OlAfgai. Section 
4 explores the evolution of the orbital population induced 
by the black hole, and Section 5 provides conclusions. 

2. MODELING TECHNIQUE 

The "adiabatic squeezing" technique for generating tri- 
axial models is discussed in detail in Paper 1. To sum- 
marize, we begin with, e.g. a spherical Hernquist (1990) 



model, having a density profile 

Ma 
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where M is the total mass and o is a scale length. We gen- 
erate an N-body realization of this spherical model using 
the multimass scheme of Sigurdsson et al. (1995), so that 
particles have a mass which is roughly inversely propor- 
tional to their pericentric radius. This technique gives a 
finer sampling of the phase space of tightly bound orbits 
than would be feasible if equal mass particles were used 
throughout. We then evolve this model while adiabati- 
cally applying a drag force first along the z axis (using the 
SCF method described below to maintain axisymmetry in 
the xy plane) and then along the y axis to yield a triax- 
ial figure. The system is rescaled after each step so that 
the scale radius of the long axis, a, is unity. The lengths 
of the other axes are correlated with the strength of the 
drag force applied in each direction. To ensure an equilib- 
rium final state, the model is evolved for several half-mass 
dynamical times without any drag forces present. 

A black hole of final mass Mbh is then grown adiabat- 
ically in one of these equilibrium models over a time fen 
according to: 



M{t) = 




for t < iBH; 

for t > tBK- 



This expression ensures smooth evolution of the black hole 
mass, i.e. that M{t) = at t = and at t = ^bh- The 
potential of the black hole is that of a softened point mass: 



$BH(r, t) 



M{t) 



BH 



where the softening parameter is set to eBH — 0.001 in the 
simulations described here. This softening sets the resolu- 
tion scale of our model and can be compared to the radius 
of influence of a 1% mass black hole, which is tbh ~ 0.03. 

To ensure that black hole growth is adiabatic, we chose 
iBH to be long compared to the orbital periods of stars in 
the inner regions of the galaxy (Sigurdsson et al. 1995). 
We varied the value chosen for fBH and determined that for 
the examples presented here a black hole growth timescale 
tBH = 20 is adequate. This can be compared with the 
dynamical times for the fiducial Hernquist sphere with 
a = 1 which are 3.14,5.96, and 8.33 at the scale radius 
(r = 1.0), effective radius (r = 1.815), and half-mass ra- 
dius (r = 2.414), respectively. After the black hole has 
reached its final mass, the model is evolved for another 20 
time units to allow the model to reach equilibrium. 

The simulations were performed with a self-consistent 
field (SCF) code (e.g. Hernquist & Ostriker 1992; Hern- 
quist et al. 1995). In this approach, the density and poten- 
tial of the galaxy are expanded in a set of basis functions, 
with the lowest order term chosen to represent the under- 
lying density profile. The expansion coefficients are deter- 
mined from the particle distribution, using n radial terms 
and {l,m) angular terms. The examples here adopted 
nmax = 10 and rumax = Imax = 6. The A^-body particles 
move in the combined potential field of the SCF expansion 
and the central black hole (Sigurdsson et al. 1997), and 
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orbital accuracy is ensured by using a high-order Hcrmitc 
integrator with variable time-stepping. The simulations 
were run on the T3E and the Blue Horizon at the San 
Diego Supercomputer Center. 

3. RESULTS OF DYNAMICAL MODELING 

We demonstrate our technique with an N=512,000 par- 
ticle Hernquist model with total mass M = 1 and initial 
axis ratios h/a = 0.85 and c/a = 0.7, measured at the half 
mass radius. To this model, wc add an adiabatically grow- 
ing black hole of final mass Mbh = 0.01. If we insist that 
this model lies on both the global and core Fundamental 
Planes (Faber et al. 1997), then the central density slope 
of 7 = 1 (neglecting the influence of the black hole) fixes 
the absolute magnitude (My ~ —21.6), effective radius 
(reff « 6 kpc), and projected central velocity dispersion 
{cFp Ri 310 km/sec) of the corresponding galaxy. Scaling 
our model to such a galaxy results in a unit length, unit ve- 
locity, and imit time corresponding to 3 kpc, 1200 km/sec, 
and 3 x 10^ years, respectively. Thus, in our model, the 
scaled half-mass dynamical timescale is 2.5 x lO'' years. 

As the black hole grows, both the cusp slope, 7, and 
projected central velocity dispersion, dp, increase. Fig- 
ure 1 shows the cusp slope (measured at ellipsoidal radius 
\ogq = —1.3) and the projected central velocity disper- 
sion of the model (measured at projected ellipsoidal ra- 
dius logQ = —2.3) as a fimction of time during and after 
black hole growth. After an initial rapid increase in 7 dur- 
ing the black hole growth phase (between < T < 20), 
the evolution in both 7 and Gp all but ceases, settling at 
the equilibrium values 7 ~ 2.05, cFp ~ 0.7. These results 
are characteristic of adiabatic black hole growth in cuspy 
galaxies and can be compared both to analytic estimates 
for adiabatic black hole growth in a spherical 7 = 1.0 
model, which predict 7 = 7/3, cTp = 0.75 (Quinlan, Hern- 
quist & Sigurdsson 1995) and to the results from N-body 
simulations of adiabatic black hole growth in a Hernquist 
sphere, where 7 ~ 2.2 and a « 0.65 (Sigurdsson, Hern- 
quist, & Quinlan 1995).^ Note that these high central 
velocity dispersions are measured on length scales smaller 
than the typical minimum spatial resolution at the dis- 
tance of nearby ellipticals; observable central velocity dis- 
persions are averaged over larger apertures and tend to 
obscure velocity cusps. 

The evolution in galaxy shape is shown in Figure 2 for 
the innermost 2%, 10%, and 50% of the galaxy (by mass). 
Prior to black hole growth, the model has a more or less 
uniform shape throughout, with a slight tendency towards 
greater triaxiality in the center (see Figure 2 at T = 
and Paper 1). As the black hole grows, the inner regions 
quickly become rounder (Figure 2, panel a); in fact, by 
the time the black hole has reached full mass, the central 
10% of the mass, corresponding to an ellipsoidal radius 
q = ^ -I- (y/by + {z/cY < 0.1 is practically spherical 
with axis ratios a : b : c = 1.0 : 0.95 : 0.92. The shape evo- 
lution in the outer regions is much less dramatic, however 
(Figure 2, panel c). Following the growth of the black hole, 
the model exhibits a marked shape gradient, becoming 
more strongly triaxial with increasing radius. This shape 
profile, essentially axisymmetric in the center and a rela- 

®The fact that our measured cusp slope is less than the analytic ■ 
finite radial range near the center, and is not the asymptotic r = \ 



tively unaltered triaxial figure further out, remains stable 
as the system settles into equilibrium from 20 < T < 40, 
even though the dynamical timescale here is much shorter 
than the integration time of the simulation. Despite the 
near sphericity at the center, this region is, perhaps sur- 
prisingly, still triaxial enough to influence the stellar or- 
bital dynamics (Statler 1987; Hunter & de Zeeuw 1992), 
at least in principle. We will return to this issue of orbits 
in the next section. 

The final state of this model features several hallmarks 
of a black hole-embedded triaxial figure. Figure 3 shows 
the properties of this object as a function of ellipsoidal ra- 
dius g at T = 40 (12.8fdyn at q = 1), well after the black 
hole has stopped growing. Figure 3a shows the 7 ~ 2 
density cusp induced by the black hole inside log g = — 1 . 
At a larger radii logg > —1, though, this plot demon- 
strates that the system retains the original Hernquist den- 
sity profile. Figure 3b shows explicitly the strong shape 
gradient in the model. Inside tbh, both the projected 
and intrinsic velocity dispersions exhibit a strong central 
cusp (panels c and d). In the outskirts, where the model 
maintains its triaxiality, the projected velocity distribu- 
tions follow ax > cTy > cTz ill accord with a triaxial model 
where a > b > c. However, inside the cusp the projected 
velocity dispersions are commensurate. Interestingly, the 
anisotropy parameter, /3 = 1 — {vf)/{v^), becomes nega- 
tive near the black hole, as many radial orbits are given a 
large tangential component. This is consistent with mod- 
els of stellar orbits around a black hole that is adiabatically 
grown, where f3 = —0.3 (Goodman & Binney 1983; Quin- 
lan et al. 1995). Exterior to the black hole's radius of 
influence, the system is radially anisotropic (/3 > 0), as 
expected for a triaxial galaxy. 

Structurally, the most dramatic change induced in the 
model by the black hole is the strong shape gradient as 
one moves inwards (from b,c = 0.9,0.8 at logg = to 
b,c = 0.95, 0.92 at logg = —2). This significant central 
roundening may arise if box, boxlet, and eccentric tube 
orbits are lost to stochasticity near the black hole (e.g. 
Norman et al. 1985; Gerhard & Binney 1985), yet sur- 
vive at larger radii. However, even the nearly spherical 
central regions may host significantly triaxial dynamics, 
as pointed out by Statler (1987) and Hunter & de Zeeuw 
(1992). Clearly, a more rigorous orbital analysis is war- 
ranted; we present this analysis in the next section. 

4. ORBITAL PROPERTIES 

In equilibrium, the structure of a galaxy and its orbital 
phase space are closely related via the Collisionless Boltz- 

mann equation. Since the intrinsic shape of an elliptical 
is dictated by the time-averaged configuration space den- 
sity of the orbits of its stars, the shape of the galaxy may 
change with time, or with radius, in response to a chang- 
ing potential. Hence information can be gleaned about 
galactic structure through the analysis of orbit families 
and their subsequent stability in our model. Both two- 
and three-dimensional analyses are useful: planar orbits 
are subject to less numerical scattering and provide sur- 
faces of section that are easy to analyze and can be com- 
pared to previous studies, while three-dimensional orbits 

of 7 = 7/3 is to be expected since our cusp slope is measured over a 
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directly trace the particle distribution and populate reso- 
nances that do not project cleanly onto a plane. 

Using the automated orbit classifier described in Paper 
1, we analyzed both the two-dimensional phase space along 
the xy and xz planes (where, as always, we take x, y, and 
z to be aligned with the major, intermediate, and minor 
axes respectively) and the three-dimensional phase space 
populated by the particles at T = 40. We remind the 
reader that the initial state of this model (the pre-black 
hole stage) is the same model analyzed in Paper 1, so a 
direct determination of the effect of the black hole can be 
made by contrasting the figures in Paper 1 to those here. 
As before, we exploited the 8-fold symmetry of the poten- 
tial to minimize noise in the particle distribution. Orbits 
were typed by the resonances in the dominant Fourier fre- 
quencies of the particle's motion in each plane (see Paper 1 
for more information on these techniques) . Following pre- 
vious work, we identified stochasticity by searching for a 
significant change in the fundamental frequency over two 
time intervals (Papc;r 1; Valluri & Merritt 1998). Each 
time interval was comprised of 50 orbital periods of a long 
axis tube at the energy of the orbit in question.^ Each 
orbit was integrated for a total of w 200 orbital times, or 
until a hard limit of T = 8000. This hard limit was chosen 
to correspond to greater than a Hubble time of evolution 
in the model when it is scaled to real galactic units for a 
galaxy at My = — 21.6. At 3i?e in such a galaxy, the total 
integration is w two Hubble times, where each integration 
interval is over 11 Gyr at this hard limit. 

We have simplified our identification of strongly chaotic 
orbits somewhat from the method used in Paper 1. Before, 
we followed the convention of Valluri & Merritt (1998) to 
define chaotic orbits. Specifically, an orbit was considered 
strongly chaotic if A/ =| /i - /2 | //qT > f„it, where /i 
and /2 are the dominant frequencies at the first and second 
time intervals, /o is the frequency of a tube about the long 
axis, T is the time between the sampled intervals, and /crit 
is the critical threshold for the onset of a strongly chaotic 
orbit. For the black hole- less model, we set the following 
chaotic threshold: /crit = 0.05-y/T/i200- The \/T time de- 
pendence was designed to describe the diffusion of chaotic 
orbits as a random-walk through phase space (Valluri & 
Merritt 1998). This threshold isolated orbits that were so 
chaotic that the orbital shape was strongly altered over 
two successive time intervals. 

In the set of calculations for Paper 1 , the precise choice 
of /crit was not important, since nearly every orbit was sta- 
ble and therefore had a negligible A/. However, adding a 
black hole amplifies the level of chaos in a model, yields a 
wide spectrum of A/s, and requires a more careful study 
of chaotic thresholds. In particular, histograms of chaotic 
subsamples preliminarily selected with the previous cri- 
terion did not show a strong (i.e. \/T) diffusion signal. 
Concerned that the rising threshold would underestimate 
the level of chaos in our models, we devised a new, sim- 
pler technique to identify chaotic orbits. In this paper, we 
identify chaos as A/ =| /i - /2 | //i > /crit, where /crit 
is a constant empirical threshold that corresponds to the 
95th percentile of the A/ distribution for the most bound 
sample (at E = —1.0) before the black hole is grown. For 



that sample, A/crit = 0.1. Choosing the threshold in this 
manner provides us with a comparative measure of the 
black hole's effect on the phase space distribution. Addi- 
tionally, the threshold was high enough to ensure that the 
orbits we have identified as chaotic did not arise from the 
miich slower diffusion of orbits due to potential noise or 
integration error. We tested the sensitivity of our results 
by sliding the threshold from 0.5/crit to 2.0/crit, and found 
that the percentage of chaotic orbits was unaffected in this 
range. In other words, chaotic orbits in this sample typ- 
ically had A/ far larger than the threshold set by any of 
these techniques. 

4.1. Planar orbit structure 

Figure 4 shows surfaces of section as a function of bind- 
ing energy for orbits confined to the xy and xz planes. As 
in Paper 1, orbits were sampled to populate the surface of 
section evenly. In the inner regions (—1.0 < E < —0.65), 
the box and boxlct phase space is entirely replaced by 
chaotic orbits. In addition, there is an easily discernible 
population of eccentric tubes that are also driven chaotic 
— on the surface of section, these orbits occupy the high 
velocity boundary at a given initial position. While the 
tubes with larger pericenters were unaffected by chaos, ec- 
centric tubes all sample more of the center of the potential, 
and are more likely to be driven chaotic. We note that the 
stable loops seen in these tightly bound slices are not a re- 
sult of motion inside the black hole's Keplerian potential; 
even at i? = — 1.0, the mean ellipsoidal radius (q) = 0.2 is 
far larger than the range of influence of the black hole. 

The outer regions of the model responded to the growth 
of the black hole in a strikingly different manner. The 
number of orbits identified as strongly chaotic fell dramat- 
ically compared to the inner regions, and these chaotic or- 
bits were randomly situated in box and boxlet phase space. 
The 2:1 resonant boxlet ("banana") is prominent and sta- 
ble in the xz plane of this model, although it does not occur 
in the original model without a black hole (sec Miralda- 
Escude & Schwarzschild 1989 or Lees & Schwarzschild 
1992 for more information on boxlet nomenclature). In the 
xy plane, the resonant orbit islands arc very weak, partic- 
ularly in the E = —0.40 slice where chaotic and regular 
resonant orbits seem to be mixed randomly throughout 
box phase space. The relative number of these scattered 
"boxlet" orbits is consistent with the orbit identification 
error of 1%. 

A common interpretation for the lack of strong chaos in 
the weakly bound orbits is that they have been integrated 
for far fewer dynamical times than the highly bound or- 
bits. For example, orbits in the innermost slice (where 
(q) = 0.2) were integrated for w 200torb, while orbits in 
the outermost slice (with (q) = 3.9) had typically reached 
only pa 50torb before the hard limit set by the analysis 
routine. Thus, under this interpretation, there has simply 
been less time to observe substantial stochastic diffusion 
in these outermost orbits, and a longer integration might 
produce the same fraction of chaotic orbits as is seen in 
the more tightly bound set. To test this hypothesis, we in- 
tegrated a subset of the lesser bound box and boxlet orbits 
{E = -0.40; (q) = 1.3) in the xz plane for « 625 orbital 
times (see Figure 5). After 200 orbital times, the orbits 

'^Hence, each integration interval corresponds to 200 dynamical times. We note that Paper 1 erroneously stated that eajch time interval 
corresponds to 12.5 dynamical times; however, this misstatement did not affect any of the conclusions in that paper. 
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were as dynamically evolved as the most tightly bound or- 
bits, but the fraction of box and boxlet phase space that 
had gone chaotic only reached 7% (versus ~ 100% of the 
box and boxlet phase space at = — 1.0). As the inte- 
gration time increased to well over two Hubble times, the 
fraction of chaotic orbits also went up. However, although 
the percentage of chaotic orbits in this subset increased 
from 2% to 16% over the course of the experiment, even 
after a scaled time of 29 Gyr, many stable centrophilic or- 
bits persist, including a population of non-resonant boxes. 
To test the extent of numerical error over long integra- 
tions, we also integrated the same subset of orbits (i.e. xz 
plane, E = —0.40) from our initial triaxial model without 
the black hole and in this experiment, only 4% of the or- 
bits were considered strongly chaotic well after two Hub- 
ble times. So, the chaos present in these weakly bound 
orbits is a real effect, not simply an artifact of the longer 
integration time. Furthermore, the long integration time 
ensured that the remaining stable box orbits were stable 
despite repeated passes through the black hole's sphere of 
influence. 

4.2. Three-dimensional orbit structure 

To study the orbital evolution induced by a central black 
hole, it is not sufficient to consider only two-dimensional 
surfaces of section. Orbits in triaxial galaxies are not sim- 
ply confined to the major and minor planes; there are 
several important 3-D resonances which serve to support 
a triaxial figure and which do not project down onto an 
identifiable resonance in a principal plane. Furthermore, 
the fact that 2-D orbits are confined to a plane dampens 
the onset of chaos. For example, in a singular logarithmic 
triaxial potential, much of the chaos present is in the 3-D 
orbits; this chaos is not reflected in the orbit families mov- 
ing in the principle planes (e.g. Papaphilippou & Laskar 
1996). 

With this in mind, we also investigated the three- 
dimensional orbital content as traced by the particle dis- 
tribution. As in Paper 1, we sorted the final particle dis- 
tribution according to binding energy and binned the dis- 
tribution into 9 slices of equal particle number. Figure 
6 shows the percentage, by mass, of tube families versus 
chaotic orbits as a function of binding energy and radius. 
Notice that the predominant orbit families in this model 
are tubes, and the precipitous drop-off in the number of 
chaotic orbits with decreasing binding energy, as seen in 
the planar sample, is also reflected here. 

Figure 7 presents the particle distribution on a frequency 
map of fy/fz versus fx/fz after 50 orbital times. Since the 
initial conditions were set by the actual N-body model, this 
frequency map is not evenly sampled. Furthermore, since 
only the initial frequencies are shown, this plot does not 
show orbital chaos. What is shown, however, are the res- 
onant regions that are populated by the particles in this 
model. In the most bound slice {E = —1.0), it is clear 
that nearly the entire population is composed of tube or- 
bits, since practically all of the 40000 orbits in this slice 
lie along the (0, 1, —1) or (1, —1, 0) resonance lines, mark- 
ing the inner long-axis tubes and short-axis tubes, respec- 
tively. In fact, the only significant area that does not lie 
along these lines also contains tube orbits; the clump near 
{fx/ fz, fy/ fz) = (0.8,0.98), contains orbits which project 
to long axis tubes in the xy-plane plane and low-order 



boxlets in another. The majority of these tube orbits are 
well outside the Keplerian potential of the black hole; these 
tubes are instead dictated by the nearly spherical stellar 
potential which extends out past logq = — 1. To under- 
score this point, we note that most of these tubes in this 
energy slice clump near the (1, 1, 1) resonance, suggesting 
that the potential is not just axisymmetric here, but nearly 
spherical (However we note that, by itself, a predominance 
of 1:1:1 resonant orbits does not necessarily mandate a 
spherical density distribution (Statler 1987; Hunter & de 
Zeeuw 1992).) 

As wo move out in the model to the least tightly bound 
particles, although the strong presence of tube orbits re- 
mains, low order boxlets are discernible, as in the (1, —2, 1) 
resonance at E = —0.65, and the (2, 0, —3) resonance at 
E = —0.2 and —0.4. These final two panels strongly re- 
semble the frequency maps for the same energy slices in 
the model without a central black hole (see Paper 1). This, 
along with the lack of shape evolution beyond logq = 0, 
indicates that the effect of the black hole does not strongly 
alter the outer regions of the model over galactic evolution- 
ary timescales. 

The general result that more loosely bound box orbits 
persist even well after a Hubble time is not new. Both an- 
alytic estimates and early numerical simulations indicated 
that while an unsoftened black hole induces stochastic- 
ity in nearly all of the box and boxlet orbits, the loss of 
triaxiality is confined to the center. For example, Ger- 
hard & Binney (1985) modeled the disruption of planar 
box orbits as a series of discrete scattering events by a 
Mbh = 0.02Mcore black hole and found that, in a Hubble 
time, box orbits extending far outside the core are unlikely 
to have experienced enough pericentric passes to have been 
substantially scattered. In a fully self-consistent N-body 
simulation of black hole growth in a post-collapse triaxial 
potential, Norman, May & van Albada (1985) showed that 
a black hole causes box phase space to be replaced with 
fully stochastic orbits on a timescalc proportional to an or- 
bit's dynamical time. This was taken to be consistent with 
the black hole scattering picture, with the caveat that the 
small number of particles in the simulation could result in 
an additional numerical source of scattering by two-body 
relaxation. 

However, recent work has suggested that the trend to- 
ward stochasticity is a more rapid and a more global 
phenomenon. Mcrritt & Quinlan (1998) conducted self- 
consistent N-body experiments of black hole growth in 
triaxial models, and determined that black hole masses 
greater than 0.003Mgai drive the galaxy toward axisymme- 
try, and that global stochasticity - where the axisymme- 
try extends well outside the half-mass radius and evolves 
on the order of an orbital period - occurred at Mbh ~ 
O.SMgai. Our more subtle and local evolution is not nec- 
essarily at odds with these results. The study conducted 
by Merritt & Quinlan involved galaxies which were more 
triaxial than ours, and which also had flat (7 = 0) cen- 
tral density profiles. The different density structure of the 
models results in different orbit populations in the two 
studies - the constant density cores of the Merritt & Quin- 
lan models are likely populated by a much higher fraction 
of box orbits than our more cuspy models (Statler 1987; 
Hunter & de Zeeuw 1992; Arnold et al. 1994). Further- 
more, the growth of a central black hole in a 7 = model 
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represents a more significant perturbation to the central 
potential than in our 7 = 1 models. Even with identical 
initial density profiles, though, a direct comparison is com- 
plex. For example, using Schwarzschild's (1979) method 
of reconstructing a 7 = 0.5 potential via a judicious selec- 
tion of orbit libraries, Valluri & Merritt (1998) found that 
slightly different initial minor axis flattenings give differ- 
ent results for the relative fraction of tightly bound and 
less bound chaotic orbits. Moreover, the degeneracy be- 
tween a given physical shape and its possible orbital con- 
tent makes any particular solution for the orbital behavior 
non-unique. It is apparent from our study, though, that 
there exists at least one set of orbit families that generate 
a triaxial figure and respond to an adiabatically growing 
central potential in such a way as to preserve its global 
shape over long timescales. 

4.3. Discussion 

It is clear that the onset of chaos in the most tightly 
bound orbit families leads to a rapid change in the inner 
structure of the model galaxies. What is less clear is the 
specific agent driving the chaos in the system — is it scat- 
tering by the central black hole itself, the effects of the 
steepened stellar density cusp, or something else entirely? 
One way to check if the stellar cusp is responsible for the 
chaos is to remove the black hole and categorize the orbits 
that result from the 7 ~ 2 frozen stellar potential. Al- 
though this is not a self-consistent galaxy model, it is a use- 
ful tool to gauge the relative importance of the black hole 
to the steep stellar cusp (and its spherical stellar poten- 
tial). We integrated a subset of particles in the E = —0.65 
energy slice in the purely stellar potential for a total of 
200torb. and classified the orbits in the manner described 
previously. In both the two- and three-dimensional orbital 
analysis, we find that rapid chaos is strongly suppressed 
when the black hole is removed. Therefore, it seems clear 
from this experiment that the steep stellar cusp is not the 
soiirce of strong chaos in these models. This result is in- 
teresting in light of the work of Merritt & Fridman (1996) 
which links the rapid onset of global chaos to the presence 
of steep (7 ~ 2) density cusps in galaxies. Earlier work by 
Schwarzschild (1993), however, demonstrated the secular 
stability of 7 ~ 2 density cusps in scale-free triaxial log- 
arithmic potentials, more in line with our results. These 
different conclusions indicate that the connection between 
ciisp slope, figure shape, and orbit stochasticity may well 
be extremely sensitive to the various orbit families initially 
present in the model. 

While the presence of a black hole is necessary for driv- 
ing chaos in our models, the mechanics of this process are 
less clear. Under a simple black hole scattering model, 
chaotic orbits should have more - or closer - pericentric 
passes within tbh than stable ones, but this is not the 
case in our sample. In fact, the degree of stochasticity de- 
pended on neither the number of pericentric passes within 
Tbh nor on the minimum Tperi- So, while the black hole 
does induce strong chaos, it seems not to do so through a 
simple scattering process. Another possibility is that chaos 
is driven by the transition between the inner spherical Ke- 
plerian potential of the black hole and the outer triaxial 
potential of the galaxy. In this picture, weakly bound box 
orbits in the outskirts of the model may spend too little 
of their orbital period in the inflection region to be signifi- 



cantly perturbed, and arc thus not driven stochastic at all. 
In our simulations, however, there was no strong correla- 
tion between A/ and the time spent within the inflection 
region, leaving open the question of what drives chaos. In 
truth, however, the graininess of the potential makes it dif- 
ficult to disentangle the subtleties of chaos from numerical 
noise in the potential, so that analytic studies are better 
suited for this question. 

In the outer regions, orbits stay regular even after re- 
peated passages near the potential center. It is possi- 
ble that many of these orbits are actually chaotic orbits 
that are "sticky" (Siopis & Kandrup 2000), with a diffu- 
sion timescale that is much longer than a few dynamical 
times. The course-grainedness of our potential seems to 
argue against this explanation; a course-grained potential 
effectively creates holes in the Arnold web (Arnold 1964) 
through which an otherwise confined orbit may escape. In 
addition, we observed the same effect in the planar orbit 
sample, and Arnold diffusion does not occur in 2-D poten- 
tials (Merritt & Fridman 1996). It seems likely, then, that 
the regularity of these orbits is real. 

Issues of orbital chaos and regularity aside, our galaxy 
model maintains its original degree of triaxiality on a 
global scale despite the presence of a massive central black 
hole. The orbit analysis described in §4 indicates that the 
model will remain relatively stable for as long as a Hub- 
ble time. It is therefore interesting to compare this stable 
galaxy model with observations of real elliptical galaxies. 

Initially, our scaling parameters were chosen to place 
the model on both the global and the core fundamental 
planes. After the evolution in structural and kinematic 
properties driven by the black hole, does the model still 
obey these relations? In the case of the global fundamen- 
tal plane, the answer is certainly yes: the changes in the 
velocity dispersion and radial density profiles occur only 
in the central regions ( logr < — 1.5, or r ^ 100 pc ), leav- 
ing the global properties of the model unchanged. For the 
core fundamental plane, the increase in central cusp slope 
( from 7=lt0 7~2) represents a significant change, 
leaving the model with a cusp slope which is perhaps a 
bit too steep for its scaled luminosity of Mb = —21.6 (see 
e.g., , Gebhardt et al. 1996). However, the 7 — Mb re- 
lationship is quite steep and shows significant scatter at 
intermediate luminosity (Gebhardt et al. 1996), so that 
this discrepancy may not be significant - ellipticals with 
Mb ^ —21 show cusp slopes of 7 = 1 — 2. Scaling our 
model to more luminous ellipticals begins to present prob- 
lems, however, since adiabatic black hole growth models 
generically predict cusp slopes steeper than that observed 
in luminous ellipticals (e.g., Bahcall & Wolf 1976; Young 
1980; Goodman & Binney 1984; Sigurdsson et al. 1995; 
Quinlan et al. 1995). 

Turning to the issue of triaxiality in ellipticals, the tri- 
axiality parameter of the model at the half-mass radius is 
T = (1 - b^)/{l - c2) = 0.5, with a flattening c/a = 0.8. 
Unfortunately, a problem arisc;s in defining the best sample 
of ellipticals with which to compare our model. To select 
for black hole embedded ellipticals, samples of radio- loud 
ellipticals might seem the best choice, but luminous ellip- 
ticals typically have a flatter cusp slope than that used 
our models. On the other hand, although a more general 
sample of ellipticals may possess a wider range of cusp 
slopes, including steeper cusps like that of our model, they 
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may also harbor black holes at a reduced rate when com- 
pared to AGN-selected samples. Bearing these; caveats in 
mind, we compare the structural properties of our model 
to the triaxiality inferred for different samples of elliptical 
galaxies. Using a combination of photometric and kine- 
matic data, Bak & Statler (2000) show that the Davies & 
Birkinshaw (1988) sample of radio galaxies is character- 
ized by a range of shapes biased towards prolate figures 
but that - as long as ellipticals are not disk-like rotators - 
triaxialities like that of our model arc common. In addi- 
tion, the typical flattening of the galaxies in that sample 
is c/a ~ 0.7, similar to, although somewhat flatter than, 
our models. Studies of larger samples of ellipticals (not 
radio selected) also show characteristic flattening similar 
to that of our model (e.g. Rydcn 1992; Trcmblay & Mcr- 
ritt 1995). While a more detailed analysis of triaxiality 
in black hole embedded ellipticals must await more com- 
plete data, at face value our model does represent well the 
observed structural properties of elliptical galaxies. 

The connection between supcrmassive black holes and 
triaxiality has important consequences for both secular 
and hierarchical galaxy evolution models. For example, 
it has been suggested that if a central black hole drives 
its host galaxy toward axisymmetry globally and rapidly, 
one possible difference between an intrinsically bright el- 
liptical (thought to be more triaxial) and a faint elliptical 
is that the stars in the faint elliptical, with their shorter 
crossing times, have had more interactions with the black 
hole and the galaxy is thus more dynamically evolved (Val- 
luri & Merritt 1998). The black hole/bulge mass relation 
can be explained in terms of galaxy evolution (Valluri & 
Merritt 1998) as well. In this scenario, spiral galaxies be- 
gin as gas-rich disks with a small triaxial bulge. The tri- 
axial potential supports fueling of the central black hole 
through material falling into it on box orbits (e.g. Nor- 
man & Silk 1983) or by gas traveling on intersecting orbits 
which drive dissipation and inflow. The black hole grows 
until a critical mass of ~ few % Afgai, whic;h breaks triaxi- 
ality and strongly curtails the gas inflow. Subsequent disk- 
disk merging can create a elliptical galaxy, and black hole 
feeding ensues in this larger triaxial bulge until the critical 
black hole mass is achieved. In both types of galaxies, the 
process is the same: once the black hole mass fraction is 
large enough to disrupt box orbits, gas inflow is sharply 
diminished. 

The black hole in our model induced axisymmetry out 
to several hundred parsecs, and resulted in a clearly ob- 
servable change in the shape and structure of the galaxy 
on these scales. Since the transformation did not take 
place globally, it is tempting to say that the black hole 
mass/bulge mass relation observed in the current galaxy 
population is not simply an artifact of gas inflow in a more 
triaxial-shaped progenitor population. However, it is not 
immediately clear how the more localized axisymmetry we 
observed would affect gas inflow and subsequent black hole 
feeding. While it is true in our globally triaxial model that 
gas inflow from outside the half mass radius would never 
be entirely cut off, the behavior of the gas once it hits the 
axisymmctric region requires detailed gas dynamical simu- 
lations. Nonetheless, it is clear that a central supcrmassive 



black hole causes dramatic and long-lasting changes in the 
host galaxy over scales well outside the region in which it 
dominates the potential. 

Finally, 21-cm observations of elliptical galaxies have 
revealed the presence of extended neutral hydrogen disks 
and rings in many ellipticals (e.g. Franx et al. 1994; van 
Gorkom & Schiminovich 1997; Hibbard et al. 2001). Sev- 
eral authors have proposed that these structures are pre- 
cursors to "disk rebuilding" in ellipticals (Schweizer 1998; 
van Gorkom 2001). If this scenario is correct, there must 
be an inward migration of gas in these systems. Triaxi- 
ality offers a mechanism to drive this inflow, but if black 
holes were to break triaxiality on large scales it is difficult 
to see how migration and disk building could occur. Our 
results alleviate these concerns; the fact that triaxiality 
is maintained in all but the inner few hundred parsecs of 
the galaxy would allow gas to move inwards on kilopar- 
sec scales to perhaps begin the process of disk formation. 
However, the modest degree of triaxiality in the model 
suggests that the rebuilding timescale may be long.^ 

5. SUMMARY 

Using numerical simulations, we have studied the 
growth of central supermassive black holes inside cuspy 
triaxial galaxies. Unlike prcndous self-cemsistent modeling 
of black hole growth in triaxial ellipticals, our calculations 
employ progenitor galaxy models which are both triaxial 
and have central density cusps typical of observed ellipti- 
cals (Holley-Bockelmann et al. 2000). Inside these progen- 
itors, we adiabatically grow a central black hole of mass 
Mbh = O.OlMgai. As the black hole grows, it induces a 
central cusp (7 ~ 2) in the stellar density profile and a 
strong roundening of the central shape of the elliptical. 
However, while the effects of the black hole do extend be- 
yond its nominal "sphere of influence" , out at an effective 
radius the galaxy figure remains largely unchanged. 

To explore the change in the orbital structure induced 
by the black hole, we use a combination of Fourier spec- 
tral classification and axis-crossing pattern recognition to 
classify the orbits present in the model. At the most 
tightly bound energies, the models are composed entirely 
of loops, short and long axis tubes, and chaotic orbits. 
These chaotic orbits comprise all of the box and boxlet 
phase space present in the original model, and even a pop- 
ulation of eccentric tubes. The outer regions predomi- 
nantly contain short axis tube orbits; however, in these 
outer regions, there still exists a modest phase space of 
boxes and boxlets which support the large-scale triaxiality 
of the system. Of these centrophilic orbit families, there 
are a substantial fraction that do not go strongly c;haotic, 
even over thousands of dynamical times. While the pres- 
ence of noise in the potential expansion limits our ability 
to detect mild chaos in the orbit populations, we believe 
that the remaining boxes and boxlets in the outer regions 
are stable enough to continue to support the global shape 
of the galaxy for a Hubble time. 

While massive black holes act as agents of chaos in the 
inner regions of our models, they provoke a more modest 
response in the outer regions than expected on the basis 

*This is particularly problematic for the very extended (i.e. many Re) rings seen in some ellipticals. In these cases, not only is the dynamical 
timescale long, but the dynamics are more likely driven by the structure of the dark matter halo than by the luminous galaxy itself (e.g., Pranx, 
van Gorkom, & de Zeeuw 1994). 
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of earlier studies (e.g. Norman et al. 1985; Merritt & 
Quinlan 1998). However, those previous studies employed 
galaxy models which had flat central density "cores" and 
were also significantly more flattened than ours, a combi- 
nation which yields a significantly higher fraction of boxes 
in the original orbital families than found in our model. 
These differences are likely the root cause of the more dra- 
matic evolution seen in previous calculations. Our galaxy 
models indicate that even moderately triaxial ellipticals 
can host central massive black holes, in agreement with ob- 
servational evidence which suggests both that black holes 
are ubiquitous (Magorrian et al. 1998) and that triaxiality 
may be common (Bak & Statler 2000). 

Our more realistic elliptical galaxy models make an ex- 
cellent tool for the simulation of several unsolved problems 



in elliptical galaxy formation and evolution, such as the 
degree of disk rebuilding by large-scale gas inflow, the ef- 
fect of satellite infall on the structure of the galaxy, the 
interaction between galaxies and their triaxial dark matter 
halos, the persistence of central cusps, and the formation 
of nuclear disks. 
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Fig. 1. — Central properties of the model galaxy during and after bla<;k hole growth. Left panel: Central density slope 7 of the model, 
measured at ellipsoidal radius logg = —1.3, as a function of time. Right panel: Central projected velocity dispersion, measured at projected 
ellipsoidal radius log Q = —2.3 as a function of time. 
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Fig. 2. — The intermediate and minor axis lengths as a function of time for particle sets binned by mass in the model galaxy. Axis lengths 
are calculated iteratively from the ellipsoidal density distribution using the moment of inertia tensor. (See Paper 1 for details.) 
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FlK. 3.- The structural and kinematic properties of the model at T = 40. Upper left: density profile. Upper right: intermediate and 
minor axis lengths as a function of ellipsoidal radius. Lower left: projected velocity dispersion along the fundamental axes, as a function of 
projected ellipsoidal radius. Lower right; true radial and tangential velocity dispersion, and velocity anisotropy parameter, as a function of 
ellipsoidal radius. 




Fig. 4. — Surfaces of section for the triaxial model at T = 40, plotted for orbital populations at different binding energies. Top: Surfaces 
of section for orbits in the xz plane. Bottom: Surfaces of section for orbits in the xy plane. Orbits are coded by color - chaotic orbits: red; 
loops: green; bananas: yellow; fish: blue; pretzels: aqua; higher resonance boxlets and pure boxes: black. This plot was created by taking an 
average of all orbit types at a particular position on the surface of section. 
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Fig. 5.— 



The percentage of chaotic orbits as a function of dynamical time for a subset of centrophilic planar xz orbits at B = 



-0.6. 
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Fig. 6. — Percentage of tubes vs. chaotic orbits in the 3-D sample as a function of binding energy. Open symbols represent tubes in general; 
open squares represent short axis tubes, open triangles show long axis tubes, and open circles show the 1:1:1 resonance, or planar loop family. 
The crosses represent strongly chaotic orbits. 
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Fig. 7. — Frequency map for the trietxial model at T = 40, plotted for orbital populations of different binding energies. The greyscale 
represents the number of orbits at a given frequency ratio. The lightest grey is 1 orbit, while black is greater than 50 orbits. The diagonal 
line corresponds to short axis tubes, and the horizontal line at fy/ fz = 1-0 corresponds to long axis tubes. 



